home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / visulztn / saoimage / saoimage.lha / histlist.c < prev    next >
C/C++ Source or Header  |  1991-06-21  |  15KB  |  429 lines
  1. #ifndef lint
  2. static char SccsId[] = "%W%  %G%";
  3. #endif
  4.  
  5. /* Module:    histlist.c (Histogram List)
  6.  * Subroutine:    make_equalized_list()    returns: void
  7.  * Copyright:    1989 Smithsonian Astrophysical Observatory
  8.  *        You may do anything you like with this file except remove
  9.  *        this copyright.  The Smithsonian Astrophysical Observatory
  10.  *        makes no representations about the suitability of this
  11.  *        software for any purpose.  It is provided "as is" without
  12.  *        express or implied warranty.
  13.  * Modified:    {0} Michael VanHilst    initial version          30 May 1989
  14.  *        {1} MVH constraints based on best history      15 Dec 1989
  15.  *        {2} JRW limit to get out of infinite loop      19 Jun 1991
  16.  *        {n} <who> -- <does what> -- <when>
  17.  */
  18.  
  19. #include <stdio.h>
  20. #include "hfiles/histeq.h"        /* define SubrangeList */
  21.  
  22. #define MAXITERATIONS 1000 /* --jrw */
  23.  
  24. #ifdef DEBUG
  25. #define HDEBUG
  26. #endif
  27.  
  28. /*
  29.  * Subroutine:    make_equalized_list
  30.  * Purpose:    Distributing levels for a subrange section of the histogram
  31.  */
  32. void make_equalized_list ( histogram, list, low_entry, high_entry,
  33.                pixel_area, color_levels )
  34.      int *histogram;
  35.      SubrangeList *list;
  36.      int low_entry, high_entry;
  37.      int pixel_area, color_levels;
  38. {
  39.   int shrink_level, stretch_level;
  40.   int max_shrink, min_stretch;
  41.   int average_area, end_area, min_area, max_area;
  42.   int levels_given;        /* levels which equalize_simply used */
  43.   int correction;        /* correction for asymptotic convergence */
  44.   int iterations;        /* avoid infinite loop --jrw */
  45.   int end_error;        /* deviation from ideal of last level */
  46.   int best_levels_over;        /* levels from an equalize worth remembering */
  47.   int best_levels_under;
  48.   int best_average_over = 0;    /* average_area used to get best_levels */
  49.   int best_average_under = 0;
  50.   static int equalize_simply();
  51. #ifdef JIGGLE
  52.   static void adjust_list();
  53. #endif
  54.  
  55.   /* else allocation distribution must first be determined */
  56.   /* run through histgram section making basic allocation and taking notes */
  57.   average_area = pixel_area / color_levels;
  58.   max_area = max_shrink = 0;
  59.   min_area = min_stretch = pixel_area + 1;
  60.   /* use 2*color_levels to allow allocation to exceed limit */
  61.   levels_given =
  62.     equalize_simply(histogram, list, 0, average_area,
  63.             low_entry, high_entry, color_levels + color_levels,
  64.             &shrink_level, &stretch_level,
  65.             &end_area, &min_area, &max_area,
  66.             &min_stretch, &max_shrink);
  67.   if( levels_given != color_levels ) {
  68.     /* make approximation for large errors */
  69.     correction = ((color_levels - levels_given) * average_area) / -100;
  70.     /* default correction if error was relatively small */
  71.     if( (correction < 2) && (correction > -2) ) {
  72.       if( color_levels > levels_given )
  73.     correction = -2;
  74.       else
  75.     correction = 2;
  76.     }
  77.     /* record what has been learned so far */
  78.     if( levels_given > color_levels ) {
  79.       best_levels_over = levels_given;
  80.       best_average_over = average_area;
  81.     } else {
  82.       best_levels_under = levels_given;
  83.       best_average_under = average_area;
  84.     }
  85.   } else
  86.     /* if first try was good enough, go with it */
  87.     correction = 0;
  88.   end_error = 0;
  89.   iterations = 0; /* --jrw */
  90.   while( (correction != 0) && (iterations++ < MAXITERATIONS) ) { /* --jrw */
  91.     average_area += correction;
  92.     levels_given =
  93.       equalize_simply(histogram, list, 0, average_area,
  94.               low_entry, high_entry, color_levels + color_levels,
  95.               &shrink_level, &stretch_level,
  96.               &end_area, &min_area, &max_area,
  97.               &min_stretch, &max_shrink);
  98.     if( levels_given != color_levels ) {
  99.       if( levels_given > color_levels ) {
  100.     /* is this the closest overshoot */
  101.     if( (best_average_over == 0) || (levels_given < best_levels_over) ) {
  102.       best_levels_over = levels_given;
  103.       best_average_over = average_area;
  104.     }
  105.     if( correction < 0 ) {
  106.       /* went too far, turn back if that is an option */
  107.       if( (correction == -1) && (end_error == 0) )
  108.         /* done if we never got to a true result */
  109.         correction = 0;
  110.       else
  111.         /* asymtotic correction to converge on solution */
  112.         correction = (correction - 1) / -2;
  113.     } else if( (best_average_under != 0) &&
  114.            (best_average_under <= (average_area + correction)) )
  115.       /* if continue in same direction, don't pass previous best */
  116.       correction = (best_average_under - average_area) - 1;
  117.       } else {
  118.     /* is this the closest undershoot */
  119.     if( (best_average_under == 0) || (levels_given < best_levels_under) ) {
  120.       best_levels_under = levels_given;
  121.       best_average_under = average_area;
  122.     }
  123.     if( correction > 0 ) {
  124.       if( (correction == 1) && (end_error == 0) )
  125.         /* done if we never got to a true result */
  126.         correction = 0;
  127.       else
  128.         /* asymtotic correction to converge on solution */
  129.         correction = (correction + 1) / -2;
  130.     } else if( (best_average_over != 0) &&
  131.            (best_average_over <= (average_area + correction)) )
  132.       /* if continue in same direction, don't pass previous best */
  133.       correction = (best_average_over - average_area) + 1;
  134.       }
  135.     } else if( (end_area > max_area) || (end_area < min_area) ) {
  136.       if( correction < -16 )
  137.     correction = -16;
  138.       else if( correction > 16 )
  139.     correction = 16;
  140.       if( end_area > average_area ) {
  141.     if( correction < 0 ) {
  142.       /* this is an overshoot */
  143.       correction = correction / -2;
  144.       if( (end_error == 0) || ((end_area - average_area) > end_error) ) {
  145.         /* go back if we might do better */
  146.         if( correction == 0 )
  147.           correction = 1;
  148.       }
  149.       end_error = end_area - average_area;
  150.     }
  151.       } else {
  152.     if( correction > 0 ) {
  153.       /* this is an overshoot */
  154.       correction = correction / -2;
  155.       if( (end_error == 0) || ((average_area - end_area) > end_error) ) {
  156.         /* go back if we might do better */
  157.         if( correction == 0 )
  158.           correction = -1;
  159.       }
  160.       end_error = average_area - end_area;
  161.     }
  162.       }
  163.     } else
  164.       correction = 0;
  165.   }
  166.   if (correction != 0) {
  167.     (void)printf("SAOimage: histeq could not converge on best distribution,");
  168.     (void)printf(" proceeding anyway\n");
  169.     (void)fflush(stdout);
  170.   }
  171.  
  172.   /* if no fit yet, force one by making last level stretch */
  173.   if( (levels_given != color_levels) && (best_average_over != 0) ) {
  174.     levels_given =
  175.       equalize_simply(histogram, list, 0, best_average_over,
  176.               low_entry, high_entry, color_levels,
  177.               &shrink_level, &stretch_level,
  178.               &end_area, &min_area, &max_area,
  179.               &min_stretch, &max_shrink);
  180.   }
  181. #ifdef JIGGLE
  182.   /* This section to directly adjust level boundaries for a better fit
  183.    * is experimental. It is not currently used, because the improvement
  184.    * is not important relative to the processing time that can occur
  185.    * in non-stable situations.  In some situations, such as when
  186.    * levels != color_levels coming out of the previous loop, a few
  187.    * passes could yield worthwhile improvements.
  188.    */
  189. #ifdef HDEBUG
  190.   (void)printf("Before adjustment: %d (between: %d %d), end: %d\n",
  191.            max_area - min_area, min_area, max_area, end_area);
  192.   (void)fflush(stdout);
  193. #endif
  194.   /* jiggle the allocations for a better fit (reduce maximum error) */
  195.   if( (end_area > max_area) && (end_area > min_stretch) ) {
  196.     /* bunch toward the end to make it smaller */
  197.     while( (end_area > max_area) && (end_area > min_stretch) ) {
  198.       average_area = min_stretch;
  199.       adjust_list(1, histogram, list, color_levels,
  200.           low_entry, high_entry, average_area,
  201.           &stretch_level, &shrink_level, &end_area,
  202.           &min_area, &max_area, &min_stretch, &max_shrink);
  203.     }
  204.   } else if( (end_area < min_area) && (end_area < max_shrink) ) {
  205.     /* bunch away from the end to make it bigger */
  206.     while( (end_area < min_area) && (end_area < max_shrink) ) {
  207.       average_area = max_shrink;
  208.       adjust_list(0, histogram, list, color_levels, low_entry, high_entry,
  209.           average_area, &stretch_level, &shrink_level, &end_area,
  210.           &min_area, &max_area, &min_stretch, &max_shrink);
  211.     }
  212.   }
  213. #ifdef HDEBUG
  214.   (void)printf("After adjustment: %d (between: %d %d), end: %d\n",
  215.            max_area - min_area, min_area, max_area, end_area);
  216.   (void)fflush(stdout);
  217. #endif
  218. #endif
  219. }
  220.  
  221. /*
  222.  * Subroutine:    equalize_simply
  223.  * Purpose:    Make a list to describe map using basic allocation method
  224.  * Note:    Allocate levels from "level" to "color_levels"
  225.  */
  226. static int equalize_simply ( histogram, histlist,
  227.                  level, average_area, low_entry, high_entry,
  228.                  color_levels, shrink_level, stretch_level,
  229.                  end_area, min_area, max_area,
  230.                  min_stretch, max_shrink )
  231.      int *histogram;
  232.      SubrangeList *histlist;
  233.      int level, average_area;
  234.      int low_entry, high_entry;
  235.      int color_levels;
  236.      int *shrink_level, *stretch_level;
  237.      int *end_area, *min_area, *max_area;
  238.      int *min_stretch, *max_shrink;
  239. {
  240.   int entry, neighbor_entry;
  241.   int area, old_area;
  242.   int err_low, err_high;
  243.   int init_next;
  244.  
  245.   /* initialize parameters */
  246.   init_next = 0;
  247.   area = 0;
  248.   histlist[level].first = low_entry;
  249.   /* reduce by one for indexing to end on last list level */
  250.   color_levels--;
  251.   /* in this pass, assign levels by simple method */
  252.   for( entry = low_entry; entry <= high_entry; entry++ ) {
  253.     if( init_next ) {
  254.       ++level;
  255.       histlist[level].first = entry;
  256.       area = 0;
  257.       init_next = 0;
  258.     }
  259.     old_area = area;
  260.     area += histogram[entry];
  261.     /* while levels last, scan till exceed average area */
  262.     if( (level < color_levels) && (area >= average_area) ) {
  263.       err_low = average_area - old_area;
  264.       err_high = area - average_area;
  265.       if( err_high < err_low ) {
  266.     /* err to excess */
  267.     histlist[level].last = entry;
  268.     histlist[level].pixel_area = area;
  269.     /* shrink option would be to exclude this entry */
  270.     histlist[level].shrink_area = old_area;
  271.     histlist[level].shrink_entry = entry - 1;
  272.     /* look for next non-zero entry to include for stretch option */
  273.     neighbor_entry = entry;
  274.     do {
  275.       ++neighbor_entry;
  276.     } while( (neighbor_entry <= high_entry) && 
  277.            (histogram[neighbor_entry] == 0) );
  278.     if( neighbor_entry > high_entry ) {
  279.       /* if at high end, exagerate to preclude this option */
  280.       histlist[level].stretch_area = 10 * area;
  281.       histlist[level].stretch_entry = high_entry;
  282.     } else {
  283.       histlist[level].stretch_area = area + histogram[neighbor_entry];
  284.       histlist[level].stretch_entry = neighbor_entry;
  285.     }
  286.       } else {
  287.     /* err short */
  288.     /* exclude this entry */
  289.     neighbor_entry = entry - 1;
  290.     histlist[level].last = neighbor_entry;
  291.     histlist[level].pixel_area = old_area;
  292.     /* including this entry is the stretch option */
  293.     histlist[level].stretch_area = area;
  294.     histlist[level].stretch_entry = entry;
  295.     /* look for previous non-zero entry for shrink option (check limit) */
  296.     while( (neighbor_entry >= low_entry) &&
  297.            (histogram[neighbor_entry] == 0) )
  298.       --neighbor_entry;
  299.     if( neighbor_entry < low_entry ) {
  300.       /* if at low end, preclude this option */
  301.       histlist[level].shrink_area = 0;
  302.       histlist[level].shrink_entry = low_entry;
  303.     } else {
  304.       histlist[level].shrink_area = old_area - histogram[neighbor_entry];
  305.       histlist[level].shrink_entry = neighbor_entry - 1;
  306.     }
  307.     /* leave this entry for the next level */
  308.     --entry;
  309.       }
  310.       /* check for new champions of excess */
  311.       if( histlist[level].pixel_area < *min_area )
  312.     *min_area = histlist[level].pixel_area;
  313.       if( histlist[level].pixel_area > *max_area )
  314.     *max_area = histlist[level].pixel_area;
  315.       if( histlist[level].stretch_area <= *min_stretch ) {
  316.     *min_stretch = histlist[level].stretch_area;
  317.     *stretch_level = level;
  318.       }
  319.       if( histlist[level].shrink_area >= *max_shrink ) {
  320.     *max_shrink = histlist[level].shrink_area;
  321.     *shrink_level = level;
  322.       }
  323.       init_next = 1;
  324.     }
  325.   }
  326.   /* mark end of list (level should be last one (color_levels - 1)) */
  327.   histlist[level].pixel_area = area;
  328.   *end_area = area;
  329.   histlist[level].last = entry - 1;
  330.   return( level + 1 );
  331. }
  332.  
  333. #ifdef JIGGLE
  334. /*
  335.  * Subroutine:    adjust_list
  336.  * Purpose:    Alter a level of the map and remake section above altered level
  337.  */
  338. static void adjust_list ( stretch_this_entry, histogram, histlist,
  339.               color_levels, low_entry, high_entry,
  340.               average_area, stretch_level, shrink_level, end_area,
  341.               min_area, max_area, min_stretch, max_shrink )
  342.      int stretch_this_entry;    /* i: stretch an entry, else shrink one */
  343.      int *histogram;
  344.      SubrangeList *histlist;
  345.      int color_levels;
  346.      int low_entry, high_entry, average_area;
  347.      int *shrink_level, *stretch_level;
  348.      int *end_area, *min_area, *max_area;
  349.      int *min_stretch, *max_shrink;
  350. {
  351.   int level, j, neighbor_entry;
  352.   int equalize_simply();
  353.  
  354.   /* make an adjustment */
  355.   if( stretch_this_entry ) {
  356.     /* stretch an entry */
  357.     level = *stretch_level;
  358.     /* shrink returns to what it was */
  359.     histlist[level].shrink_area = histlist[level].pixel_area;
  360.     histlist[level].shrink_entry = histlist[level].last;
  361.     /* current was stretch option before */
  362.     histlist[level].pixel_area = histlist[level].stretch_area;
  363.     histlist[level].last = histlist[level].stretch_entry;
  364.     /* find the next non-zero entry for the new stretch option */
  365.     neighbor_entry = histlist[level].last;
  366.     do {
  367.       ++neighbor_entry;
  368.     } while( (neighbor_entry <= high_entry) &&
  369.          (histogram[neighbor_entry] == 0) );
  370.     if( neighbor_entry > high_entry ) {
  371.       /* don't go off the edge */
  372.       histlist[level].stretch_area *= 10;
  373.       histlist[level].stretch_entry = high_entry;
  374.     } else {
  375.       histlist[level].stretch_area += histogram[neighbor_entry];
  376.       histlist[level].stretch_entry = neighbor_entry;
  377.     }
  378.   } else {
  379.     /* shrink an entry */
  380.     level = *shrink_level;
  381.     /* stretch option becomes what it was */
  382.     histlist[level].stretch_area = histlist[level].pixel_area;
  383.     histlist[level].stretch_entry = histlist[level].last;
  384.     /* current was shrink option before */
  385.     histlist[level].last = histlist[level].shrink_entry;
  386.     histlist[level].pixel_area = histlist[level].shrink_area;
  387.     /* look for last non-zero entry to exclude for new shrink option */
  388.     neighbor_entry = histlist[level].last;
  389.     while( (neighbor_entry >= low_entry) &&
  390.        (histogram[neighbor_entry] == 0) )
  391.       --neighbor_entry;
  392.     if( neighbor_entry < low_entry ) {
  393.       /* don't go off the edge */
  394.       histlist[level].shrink_area = 0;
  395.       histlist[level].shrink_entry = low_entry;
  396.     } else {
  397.       histlist[level].shrink_area -= histogram[neighbor_entry];
  398.       histlist[level].shrink_entry = neighbor_entry - 1;
  399.     }
  400.   }
  401.   /* remake list with this change */
  402.   *min_area = *min_stretch = average_area * 3;
  403.   *max_shrink = *max_area = 0;
  404.   /* retake notes up to point of change */
  405.   for( j = 0; j <= level; j++ ) {
  406.     if( histlist[j].pixel_area < *min_area )
  407.       *min_area = histlist[j].pixel_area;
  408.     if( histlist[j].pixel_area > *max_area )
  409.       *max_area = histlist[j].pixel_area;
  410.     if( histlist[j].stretch_area <= *min_stretch ) {
  411.       *min_stretch = histlist[j].stretch_area;
  412.       *stretch_level = j;
  413.     }
  414.     if( histlist[j].shrink_area >= *max_shrink ) {
  415.       *max_shrink = histlist[j].shrink_area;
  416.       *shrink_level = j;
  417.     }
  418.   }
  419.   if( (level < color_levels) && (histlist[level].last < high_entry) )
  420.     /* make rest of list in normal way */
  421.     (void)equalize_simply(histogram, histlist, j, average_area,
  422.               histlist[level].last + 1, high_entry, color_levels,
  423.               shrink_level, stretch_level,
  424.               end_area, min_area, max_area,
  425.               min_stretch, max_shrink);
  426. }
  427. #endif
  428.  
  429.